
scaledMoments <- function(moments) {
	m1 <- moments[1]
	m2 <- moments[2]
	m3 <- moments[3]
	m4 <- moments[4]
	var_z <- -m1^2 + m2
	sd_z <- sqrt(var_z)
	skew_z <- (m3 - 3*m1*sd_z^2 - m1^3)/(sd_z^3)
	kurt_z <- (-3*m1^4 + 6*m1^2*m2 - 4*m1*m3 + m4)/(sd_z^4)
	return(c(m1, var_z, skew_z, kurt_z))
}

set.seed(2138)

n <- 100000
mean_z <- 10
sd_z <- 3
S <- sd_z
Z <- rnorm(n, mean_z, sd_z)
X <- Z + rnorm(n, 0, sd_z)
b0 <- 10
b1 <- 3
e <- c(rnorm(n/2, 0, sd = 1),  rnorm(n/2, 0, sd = 6))
W <- 10 + b1*Z + e
Y <- W + rnorm(n, 0, sd_z)

private_model <- lm(W ~ Z)

data <- data.frame(rbind(c(sd_z, sd_z), cbind(X, Y)))
dp_model <- lmdp(Y ~ X, data = data)

beta <- dp_model$beta_tilde
dp_res <- Y - cbind(1, X)%*%beta
err_X <- c(0, sd_z^2)


plot <- ggplot() + 
  geom_density(aes(x = dp_res, fill= 'Observed residuals'), alpha = .7) + 
  geom_density(aes(x = e, fill = 'True errors'), alpha = .5) +
   theme_bw() + 
   theme(legend.position = c(0.9, 0.9), 
       legend.justification =  c(0.9, 0.9),
       panel.grid.major = element_blank(), 
       panel.grid.minor = element_blank()) + 
  labs(x = '', fill = '') + scale_fill_manual(values = c('darkblue', 'orange'))
  

ggsave(plot, file = 'plots/figure5.pdf', width = 6, height = 4)

# Observed moments
obs_moments <- c(moments::skewness(dp_res), moments::kurtosis(dp_res))

# Estimated moments
est_moments <- scaledMoments(Rmoments(X = as.numeric(dp_res), S = sqrt(beta[2]^2*S^2 + S^2), R = 4, n = n)[, 1])[3:4]

table2 <- as.data.frame(cbind(obs_moments, est_moments))

write.csv(table2, 'table3.csv')



